The three body problem solution using Runge-Kutta method 

Juan Urrutia-Mustapha Bousackla, Universitat Autonoma de Bercelona 

January 18, 2019 


Abstract 

In order to solve the 3 body problem, it’s necessary to compute the movement via a numerical 
method such as Runge-Kutta. Using these method, it’s shown the movement of the the Earth, Sun and 
an asteroid. 


1 Introduction 

The three body problem has been an important challenge in physics history. The first one to come up with 
a possible solution was Leonhard Euler who found three families of periodic solutions as long as the three 
masses were on the same line. About twenty years later on 1772, Lagrange found a family of stable orbits 
for the three masses, made of ellipses. As response to the 1884 competition that Oscar II of Sweden made 
for solving the problem, Henri Poincare answered with the recurrence theorem which made him won. The 
theorem states that for certain systems after a big but finite time they advance to the initial state. It was not 
until 1912 that Sundman found a convergent solution for systems that have a non zero angular momentum. 
But the solution is not for practical use. The most used method today is a complex numerical method. 


2 Numerical method 


2.1 Normalization 


In order to avoid the truncation error and use more of the computer’s memory the variables are normalized 
using the third Kepler law applied to the Earth Sun system: 


And the Newton’s gravity law: 
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By combining both equations it’s obtained: 
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Using the following normalization: 
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2.2 Runge-Kutta 

The use of simpler method as Euler will give too much error. The Runge-kutta method is used for solving 
differential equations as ec(3). It’s needed to transforms those ecuations in something like this: 

y' = f{t,y ) 
y(t 0 ) = y a 

In every h time step we define the k\ , k 2 , &3 and k± as follows. 

ki = f{t,y) 

k 2 = f(t+ ^,y + ki^) 

, „/ h , h . 

k3 = f(t+ 2 ^y + k 2 ~) 

kA = f(t+ 7^,y + k 2 h ) 


So to compute the solution one step forward, we define that: 

Vn +1 = Vn + tt(&i + 2 fc 2 + 2 /C 3 + k 4 ) 
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2.3 Runge-Kutta applied to 3 body problem 

The formulation for the three body problem applies in the following way: 

v t+\ = v\ + — (k[ + 2 fc 2 + 2fcg + k l 4 ) 
r'i+i = < + f {K\ + 2 \K\ + 2 K\ + K\) 
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For every time step it’s necessary to compute every K and k in the following way. Where the sub index ”i” 
refers to the component on the 3 dimensional space. 
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-)> 
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a l( r i + v fdt) = K l 2 
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Consult the complete code on the annex. 

3 Earth Sun and asteroid system 

3.1 Data 

Using the normalized values for the system it’s obtained: 

dearth = 3 * 10 ' 6 

J — I 
'earth x 

And the following data for the asteroid: 


^max 11-59 

^min 0-6 
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3.2 Results and precision 

The accumulative error of the 4 orden Runga-Kutta method is of 0(dt 4 ). In the system earth sun asteroid 
it’s obtained and error of about 0.0005 in one year(Fig.2). This is the result obtained for an arbitrary asteroid 



Figure 1: The result for an arbitrary asteroid 

with the specified characteristics. We can see the movement of the earth along the sun in one year and the 
eccentric orbit of the asteroid. 

3.3 Conclusions 

Runge-Kutta is a good method for solving newton’s gravity. By changing the original data of the program 
developed it’s possible to solve many problems such as: the moon-earth-sun system, the Kepler asteroid 
movement or the Oumuamua approximation to earth. Due to the needed of calculating many parameters in 
every time iteration it’s difficult but not impossible to generalize the 3 body problem to the n-body problem 
via this method. 
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4 Code 


#i n c1u d e <stdio.h> 

#i n c 1 u d e <math. h> 

#i n c1u d e <stdlib.h> 

#in c1u d e < s t rin g . h> 

#include <fstream> 

#include <iostream> 

//Where n is the numbre of iterations 

^define n 10000 

//These are the mases 

^define ml 1 

#define m2 0.000003 

#define m3 0 

//Time is normalized to 1 year 
^define t 1.0 
int main () { 

std : : ofstream fout ; 

fout.open (”tres_cuerpos . txt”) ; 

double G; 

int i ; 

int j ; 

G=4*(3.14159265) *(3.14159265) ; 

//In here we define the initial conditions 

double dt=t/n; 

double rl [ 3 ] [ n + 1]; 

double r2 [ 3 ] [ n +1]; 

double r3 [ 3 ] [ n + 1]; 

double vl [ 3] [ n+1]; 

double v2 [ 3 ] [ n+1]; 

double v3 [ 3 ] [ n+1]; 

double clld2[n + l]; 

double dld3[n + l]; 

double d2d3[n + l]; 

double kl [ 3 ] [ 3 ] ; 

double k2[3][3 ] ; 

double k3 [ 3 ] [ 3 ] ; 

double k4 [ 3 ] [ 3 ] ; 

double Kl [ 3 ] [ 3 ] ; 

double K2 [ 3 ] [ 3 ] ; 

double K3 [ 3 ] [ 3 ] ; 

double K4 [ 3 ] [ 3 ] ; 

//The position of the first body 
rl[0][0]= 0; 
rl [1] [0] = 0; 
rl [2] [0] = 0; 

//The position of the second body 
r2[0][0] = 1; 
r2[l][0]=0; 
r2 [2] [0] =0 ; 

//The position of the third body 
r3 [0] [0] = 1.00026; 
r3 [1] [0] = 0; 
r3 [2] [0] =0; 

//Velocity of the first body 



vl[0][0]=0; 

vl[l][0]=0; 

vl[2][0]=0; 

//Velocity of the second body 
v2 [0] [0] =0; 

v2[l][0]=2*(3.141516) ; 
v2 [2] [0] =0; 

//Velocity of the third body 
v3 [0] [0] =0; 
v3 [1] [0] =0; 
v3 [2] [0] =0; 

//These are the titles for the file 
fout « std : : fixed ; 

fout « ” Three body problem” «std::endl; 
fout « ” ” «std : : endl ; 

//Column names 

fout « ” Position m2(x) Position m2(y) ” «std : : endl ; 

//Temporal iteration 
for (i=0; i<n; i++) { 

/ /What 

fout «r2 [0] [ i]<<” ”«r2 [ 1 ] [ i]<<std : : endl; 

//These is the K for the first body 
//Kl 

dld2 [ i ]=pow(pow( r 1 [0] [ i] — r2 [0] [ i] , 2)+pow( rl [1] [ i] — r2 [1] [ i] , 2)+pow( rl 

[2] [ i]-r2 [2] [ i] ,2) ,3/2); 

dld3 [ i ]=pow(pow( rl [0][i] — r3[0][i] , 2)+pow( rl [1] [ i] — r3 [1] [ i] , 2)+pow( rl 
[2] [ i] r3 [2][ i] , 2) ,3/2) ; 
for (j =0; j<=2; j++) { 

Kl [ j ] [0] =G*m2/dld2[i]*(r2 [j ] [ i ] — r 1 [j ] [ i] )4G*m3/ dld3 [ i]*(r3 [j ] [ i] — 
rl [j ] [ i]) ; 

kl[j][0] = vl [ j ] [ i]; 

} 

//K2 

dld2 [ i ]=pow(pow( rl [0][i] + kl[0][0]*dt /2—r2 [ 0 ] [ i ] , 2)+pow( rl [ 1 ] [ i ] + kl 
[1] [0] * dt/2-r2 [1] [ i ] , 2)+pow(rl [2] [ i] + kl [2] [0] * dt/2—r2 [2] [ i ] ,2) 
,3/2); 

dld3 [ i ]=pow(pow( rl [0] [ i] + kl [0] [0]* dt/2—r 3 [0] [ i ] , 2)+pow( rl [1] [ i] + kl 
[1] [0] * dt/2-r3 [1] [ i ] , 2)+pow(rl [2] [ i] + kl [2] [0] * dt/2—r3 [2 ] [ i ] , 2) 
,3/2); 

for (j =0; j<=2; j++) { 

k2 [ j ] [0] = kl [j ] [0] +K1 [ j ] [0]* dt/2; 

K2 [ j j [0] =G*m2/dld2 [ i ] * ( r2 [ j ] [ i] — kl [ j ] [0 ] * dt/2—rl [ j ] [ i ]) 4G*m3/dld3 [ 
i ] * ( r 3 [ j ] [ i ] — kl [ j ] [0] * dt/2—rl [ j ] [ i ]) ; 

} 

//K3 

dld2 [ i ]=pow(pow( rl [0] [ i] + k2 [0] [0]* dt/2—r 2 [0] [ i ] , 2)+pow( rl [1] [ i] + k2 
[1] [0] * dt/2-r2 [1] [ i ] , 2)+pow(rl [2] [ i] + k2 [2] [0] * dt/2-r2 [2 ] [ i ] ,2) 
,3/2); 

dld3 [ i ]=pow(pow( rl [0] [ i] + k2 [0] [0]* dt/2—r 3 [0] [ i ] , 2)+pow( rl [1] [ i] + k2 
[1][0]* dt/2—r3 [1][i] , 2)+pow(rl [2][i] + k2 [2][0]* dt/2—r3 [2][i] , 2) 
,3/2); 

for (j =0; j<=2; j++) { 

k3 [ j ] [0] = k2 [ j ] [0] +K2 [ j ] [ 0 ] * dt / 2; 

K3 [ j ] [0] = G*m2/dld2 [ i ] * ( r2 [ j ] [ i] — k2 [ j ] [0 ] * dt/2—rl [ j ] [ i ]) 4G*m3/dld3 [ 
i ] * ( r 3 [ j ] [ i] —k2 [ j ] [0] * dt/2—rl [ j ] [ i ]) ; 
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} 

//K4 

dld2 [ i ]=pow(pow( r 1 [ 0 ] [ i ] + k3 [0] [0] * dt—r2 [ 0 ] [ i ] , 2)+pow( rl [ 1 ] [ i ] + k3 

[1][0] * dt—r2 [1] [i] , 2)+pow(rl [2][i] + k3[2][0] *dt-r2[2][i] ,2) ,3/2); 
dld3 [ i ]=pow(pow( rl [ 0 ] [ i ] + k3 [0 ] [0] * dt—r3 [ 0 ] [ i ] , 2)+pow( rl [ 1 ] [ i ] + k3 

[1] [0]* dt—r3[1][i] , 2)+pow(rl [2][i] + k3[2][0] *dt-r3[2][i] , 2),3/2); 
for (j =0; j<=2; j++) { 

k4[j][0] = k3[j][0]+K3[j][0]*dt; 

K4 [ j ] [0] =G*m2/dld2[i]*(r2 [ j ] [ i] — k3 [ j ] [0] * dt—rl [ j ] [ i ])-fG*m3/dld3 [ i 
]*( r 3[j][i] — k3[j][0]* dt-rl [j ][i]) ; 

//These is to compute the position and the speed 
vl[j ] [ i+l]=vl [ j ] [i]+0.167*(Kl[j][0] + 2*K2[j][0] + 2*K3[j][0]+K4[j 
] [0]) *dt; 

rl [j ][i+l] = rl [j ] [i]+0.167*(kl[j][0] + 2*k2[j][0] + 2*k3[j][0] + k4[j 
] [0]) *dt; 

} 

//K for the second body 

//Kl 

dld2 [ i ]=pow(pow( rl [0][i] — r2[0][i] , 2)+pow( rl [l][i] — r2[l][i] , 2)+pow( rl 

[2] [i] r2 [2][i] ,2) ,3/2); 

d2d3[i ]=pow(pow( r2 [0][i ] — r 3 [0][i] , 2)+pow(r2 [l][i] — r3[l][i] , 2)+pow(r2 
[2][i] — r3 [ 2][i] , 2) ,3/2) ; 
for (j =0; j<=2; j++) { 

Kl [ j ] [1] =G*m3/d2d3 [ i ] * ( r3 [ j ] [ i ] — r 2 [ j ] [ i ]) 4G*ml/dld2 [i]*(rl[j][i] — 

r - [ .i ] [ i!); 

kl [ j ] [ 1] = v2[j ][ i ] ; 

} 

//K2 

dld2 [ i ]=pow(pow( rl [0] [ i ] — r 2 [0] [ i ] — kl [0] [1]* dt/2, 2)+pow( rl [ 1 ] [ i ] — r 2 
[1] [ i] — kl [1] [l]*dt/2, 2)+pow( rl [ 2 ] [ i ] — r 2 [ 2 ] [ i ] — kl [ 2 ] [ 1 ] =t= dt / 2 ,2) 

,3/2); .. 

d2d3 [ i ]=pow(pow( r2 [0][i] + kl[0][l]*dt /2—r3 [ 0 ] [ i ] , 2)+pow( r2 [ 1 ] [ i ] + kl 
[1] [1] * dt/2-r3 [1] [ i ] , 2)+pow(r2 [2] [ i] + kl [2] [ 1 ] * dt/2-r3 [2 ] [ i ] , 2) 
,3/2); 

for (j =0; j<=2; j++) { 

k2[j][l] = kl[j][l] +K1 [ j ] [ 1 ] * dt / 2; 

K2 [ j j [ 1 ] = G*m3/d2d3 [ i ] * ( r3 [ j ] [ i] — r2 [ j ] [ i] — kl [ j ] [ 1 ] * dt /2)4G*ml/dld2 [ 
i]*(rl [j ] [ i]-r2 [j ] [ i] —kl [ j ] [ 1 ] * dt/2) ; 

} 

//K3 

dld2 [ i ]=pow(pow( rl [0] [ i ] — r 2 [0] [ i ] — k2 [0] [1]* dt/2, 2)+pow( r 1 [ 1 ] [ i ] — r 2 
[1][i] — k2[1][1]* dt/2, 2)+pow(rl [2] [ i] —r2 [2][i] —k2[2][1]* dt/2,2) 
,3/2); 

d2d3 [ i]=pow(pow(r2 [0] [ i] + k2 [0] [1]* dt/2—r 3 [0] [ i ] , 2)+pow( r2 [1] [ i] + k2 
[1] [1] * dt/2-r3 [1] [ i ] , 2)+pow(r2 [2] [ i] + k2 [2] [ 1 ] * dt/2-r3 [2 ] [ i ] , 2) 
,3/2); 

for (j =0; j<=2; j++) { 

k3 [ j ] [ 1] = k2[j][1] + K2[j ][1]* dt /2; 

K3 [ j ] [ 1] = G*m3/d2d3 [ i ] * ( r3 [ j ] [ i ] — r2 [ j ] [ i ] — k2 [ j ] [ 1 ] * dt /2)4G*ml/dld2 [ 
i ] * ( r 1 [ j ] [ i ] — r 2 [ j ] [ i] — k2 [ j ] [1] * dt/2) ; 

! 

dld2[i ]=pow(pow( rl [ 0][i] — r2 [0][i] — k3[0][1]* dt , 2)+pow(rl [1][i] — r2 [1] [ i 
] —k3[1][1] * dt, 2) +pow(r1 [2][i] — r2[2][i] — k3[2][1 ] * dt ,3/2); 
d2d3 [ i ]=pow(pow( r2 [ 0 ] [ i ] + k3 [0] [ 1 ] * dt—r3 [ 0] [ i ] , 2)+pow( r2 [ 1 ] [ i ] + k3 

[1][1]* dt—r3[1][i] , 2) +pow(r 2 [2][i] + k3[2][1] *dt-r3[2][i] , 2),3/2); 
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for (j =0; j<=2; j++) { 

//K4 

k4 [ j ] [!] = k3 [ j ] [1]+K3[ j ] [1] * dt ; 

K4 [ j ] [ 1 ] = G*m3/d2d3[i]*(r3 [ j ] [ i ] — r2 [ j ] [ i ] — k3 [ j ] [ 1 ] * dt )-tG*ml/dld2 [ i 
]*(rl [j ] [ i ] — r 2 [j ] [ i] — k3 [ j ] [1] * dt) ; 

} 

//Compute the second body position and speed 
for (j =0; j<=2; j++) { 

v2[j ] [ i +I]=v2 [ j ] [i]+0.167*(Kl[j][l] + 2*K2[j][l]+2*K3[j][l]+K4[j 
] [1]) *dt; 

r 2 [ j ] [ i +I] = r2 [ j ] [ i ]+0.167* (kl [ j ] [1]+ 2* k2 [ j ] [1]+ 2* k3 [ j ] [1] + k4 [ j 
] [1]) *dt; 

} 

//Compute the k for the third body 

//Kl 

dld3 [ i ]=pow(pow( rl [0] [ i] — r3 [0] [ i] , 2)+pow( rl [l][i] — r3[l][i] , 2)+pow( rl 
[2][i] — r3 [ 2][i] , 2) ,3/2); ~ 

d2d3 [ i ]=pow(pow( r2 [0] [ i] — r3 [0] [ i] , 2)+pow( r2 [ 1 ] [ i] — r3 [ 1 ] [ i ] , 2)+pow( r2 
[2][ i] r3 [2][ i] , 2) ,3/2) ; 
for (j =0; j<=2; j++) { 

Kl [ j ] [2] = G*ml/dld3 [i]*(rl[j][i] — r3[j][i] )-fC*m2/d2d3[i]*(r2 [j ] [ i] — 
r 3 [ j ] [ i ]) ; 
kl [ j ] [2]= v3 [ j ] [ i ] ; 

} 

//K2 

d2d3 [ i ]=pow(pow( r2 [0] [ i ] — r 3 [0] [ i] — kl [0] [2]* dt/2, 2)+pow( r2 [ 1 ] [ i ] — r 3 
[1] [ i] — kl [1] [2] * dt/2, 2)+pow(r2 [2] [ i]-r3 [2] [ i]-kl [2] [2] * dt/2, 2) 
,3/2); 

dld3[i]=pow(pow(rl [0][i] — r3 [0][i] — kl[0][2]* dt/2, 2)+pow(rl [ 1 ][i] — r3 
[1][i ] — kl[1][2]* dt/2, 2)+pow(rl [2] [ i]-r3 [2][i] kl [2][2]* dt/2 , 2) 
,3/2); 

for (j =0; j<=2; j++) { 

k2[j][2] = kl[j][2]+Kl[j][2]*dt/2; 

K2 [ j ] [2] = G*ml/dld3 [ i ] * ( rl [ j ] [ i] — r3 [ j ] [ i] — kl [ j ] [2 ] * dt /2)-fG*m2/d2d3 [ 
i]*( r 2 [ j ][i] — r 3 [ j ] [ i] —kl [ j ] [2 ] *dt/2); 

} 

//K3 

d2d3 [ i ]=pow(pow( r2 [0] [ i ] — r 3 [0] [ i ] — k2 [0] [2]* dt/2, 2)+pow( r2 [ 1 ] [ i ] — r 3 
[1][i ] — k2[1][2]* dt/2, 2)+pow(r2 [2][i] — r 3 [2][i]-k2 [2][2]* dt/2 , 2) 

•3 2): 

dld3 [ i ]=pow(pow( rl [0] [ i ] — r 3 [0] [ i ] — k2 [0] [2]* dt/2, 2)+pow( r 1 [ 1 ] [ i ] — r 3 
[1] [ i ] — k2 [1] [2] * dt/2, 2)+pow(rl [2] [ i]-r3 [2] [ i]-k2 [2] [2] * dt/2 , 2) 
,3/2); 

for (j =0; j<=2; j++) { 

k3[j][2] = k2[j][2]+K2[j][2]*dt/2; 

K3 [ j j [2] = G*ml/dld3 [i]*(rl [ j ] [ i] — r3 [ j ] [ i] — k2 [ j ] [2 ] * dt /2)4G*m2/d2d3 [ 
i ] * ( r 2 [ j ] [ i ] — r 3 [ j ] [ i]-k2 [ j ] [2] * dt/2) ; 

} 

//K4 

d2d3 [ i ]=pow(pow( r2 [ 0 ] [ i] — r3 [ 0 ] [ i] — k3 [0] [ 2 ] * dt , 2)+pow( r2 [1] [ i ] — r3 [1] [ i 
]-k3 [1] [2]* dt/2, 2) +pow ( r 2 [2] [ i]-r3 [2] [ i]-k3 [2] [2] * dt/2, 2),3/2); 
dld3 [ i]=pow(pow( rl [0] [ i] — r3 [0] [ i] — k3 [0] [2] * dt , 2)+pow( rl [1] [ i] — r3 [1] [ i 
]-k3 [1] [2]* dt/2, 2) +pow ( r 1 [ 2 ] [ i]-r3 [2 ] [ i]-k3 [2] [2] * dt/2 , 2),3/2); 
for (j =0; j<=2; j++) { 

k4[j][2] = k3[j][2]+K3[j][2]*dt; 
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K4 [ j ] [2] =G*ml/dld3[i]*(rl [j ] [ i ] — r 3 [j ] [ i ] — k3 [ j ] [ 2 ] * dt )+G*m2/d2d3 [ 
]*(i'2 [ j ] [ i]-r3 [ j ] [ i] —k3 [ j ] [2] * dt) ; 


for (j =0; j<=2; j++) { 

//Calculate the speed and position for the third body 
v 3 [ j ] [ i+l]=v3[j ] [i]+0.167*(Kl[j][2]+2*K2[j][2] + 2*K3[j][2]+K4[j 
] [2]) *dt; 

r 3 [ j ] [ i+l] = r3 [ j ] [ i ]+0.167* (kl [ j ] [2]+2* k2 [ j ] [2]+2* k3 [ j ] [2] + k4 [ j 
] [2]) *dt; 
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